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Abstract. We study the effects of temperature and sliding velocity on 
superlubricity in numerical simulations of the Frenkel-Kontorova model. We show 
that resonant excitations of the phonons in an incommensurate sliding body lead 
to an effective friction and to thermal equilibrium with energy distributed over 
the internal degrees of freedom. For finite temperature, the effective friction can 
be described well by a viscous damping force, with a damping coefficient that 
emerges naturally from the microscopic dynamics. This damping coefficient is 
a non-monotonic function of the sliding velocity which peaks around resonant 
velocities and increases with temperature. At low velocities, it remains finite and 
nonzero, indicating the preservation of superlubricity in the zero-velocity limit. 
Finally, we propose experimental systems in which our results could be verified. 
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1. Introduction 

Recent developments in surface scanning probes have stimulated a large interest in 
understanding and using the phenomenon of friction at the nanoscale [U [2] . One of 
the most fascinating aspects of this field is the possibility of superlubricity, namely 
the possibility to slide two surfaces onto each other with vanishingly small friction. 
Several experiments report such a low- friction state [3l IH [5] . The crucial ingredient 
of this phenomenon is the geometrical incommensurability of the sliding crystal faces. 
Other aspects like roughness of the sliding surfaces |6j or electronic excitations during 
sliding U are not discussed here. The key model to study incommensurate systems 
is the Frenkel-Kontorova model (FK-model) [51 (shown in the top panel of figure [T]), 
which describes a chain of atoms with period a on a potential profile with period b. 
The case where the ratio b/a is not a rational number identifies an incommensurate 
system with no finite global periodicity. 

Aubry has shown that a structural transition occurs at a critical coupling Ac of the 
chain to the substrate potential [51 [TO]- Below this critical coupling, the static friction 
vanishes because, while moving, the atoms of the chain always explore all possible 
positions on the substrate and thereby average out the effective corrugation and 
resulting friction force. One of the empirical laws that describes macroscopic friction 
states that kinetic friction is lower than static friction [T . Therefore it is important to 
establish whether vanishing static friction indeed implies vanishing kinetic friction also 
at the nanoscale. In the paper where they coined the term superlubricity, Shinjo and 
Hirano [11] studied the kinetic contribution to friction in the FK model and proposed 
a phase diagram (bottom panel figured]) with a frictional and a superlubric region as 
a function of the coupling A and the velocity of the centre of mass of the chain. In 
the superlubric regime, once in motion, the chain would slide indefinitely without any 
kinetic friction but with a recurrent exchange of kinetic energy between the centre of 
mass (CM) and a single internal mode related to the periodic modulation due to the 
substrate. According to their findings, by increasing the velocity of the chain for weak 
coupling (A < Ac, no static friction), namely along a vertical line in the phase diagram 
of figured] the system would first acquire friction and then go back to a superlubric 
state at higher velocities. The proposed behaviour implies a surprising sudden change 
to a frictional state as soon as the system is set into motion. 

This result has later been criticized [12] by showing that the initial superlubric 
behaviour is destroyed with time, since the recurrent oscillation of the CM gives rise 
to the resonant parametric excitation of acoustical, long wavelength, vibrations in the 
chain. This phenomenon leads to a high friction regime where all the initial kinetic 
energy of the CM is eventually converted into internal motion that can be interpreted 
as heat. The resonant excitation of acoustical modes is however weaker at lower 
velocities and goes naturally towards the vanishing friction of the static limit. This 
result agrees well with the fact that most experimental observations of quasi-vanishing 
friction have been at the extremely low velocities of atomic force microscopes (AFM), 
of the order of /xm/s. Indeed, the term superlubricity is currently used to indicate in 
general a low friction sliding rather than in the original meaning of reference |11| . 

Here we study how superlubricity is affected by temperature and velocity. First we 
show that the dynamical mechanism identified in reference [12j leads to a system in 
thermal equilibrium that satisfies energy equipartition over the internal vibrations, 
yielding a well defined temperature. We then show that at finite temperature 
the system displays a damped dynamics that allows to calculate the damping 
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coefficient from tlie microscopic description rather than by assuming its presence 
phenomenologically. We find that this damping is not a monotonic function of the 
sHding velocity, as it peaks close to specific phonon resonances. Away from these 
resonances and particularly in the limit of low velocity, the effective damping is finite, 
and thus at low velocities the friction vanishes, preserving superlubricity. 

The paper is organized as follows. In section [2] we describe the model and its 
dynamics, in section |3] we review how kinetic friction arises and show how the system 
evolves into thermal equilibrium. In section |4] we show that, at thermal equilibrium, 
friction can be described as a viscous damping that can be calculated from the model. 
We also compare our results to damped driven models [13]. In section [5] we suggest 
experimental systems that could confirm the results presented in this paper and in 
section [S] we summarize our findings. 

2. Model: geometry and dynamics 

The FK- model [SJ consists of a harmonic linear chain of particles of mass m with period 
ao, subjected to an external periodic potential with period Os and amplitude Uq/^tt. 
The spring constant between the particles is K. The ratio of the depth of the potential 
to the spring constant of the chain defines a coupling parameter: A = Uo/{Ka^). In 
the context of friction the particles represent the atoms of a sliding body and the 
potential represents a rigid solid substrate. The model is depicted schematically in 
the top panel of figure [TJ 

If the ratio of the parameters ao/og is an irrational number, the system is called 
incommensurate, otherwise the system is called commensurate. The incommensurate 
case is the most interesting, because for two arbitrary surfaces in contact a common 
periodicity is not likely. 

The ground state is determined by the competition of the harmonic interaction 
between the particles and the interaction of the particles with the periodic substrate, 
which is controlled by A. For incommensurate systems there exists a structural phase 
transition [9l [10]: for A < Ac the system is in a so called floating phase and there 
is no static friction. If A > Ac the system is in the pinned phase and exhibits static 
friction. The value of Ac is largest when the ratio between the length scales is equal 
to the golden mean, Tg, ao/as = Tg = (1 + y5)/2 = 1.618033989 . . . [M], for which 
Ac = 0.1546... jl5j . In the numerical implementation, one may approximate the 
golden mean by a ratio of successive Fibonacci numbers, Tg — lim„_>.oo ^n+i/^n. 
Here, all calculations are done for n = 12, i.e. N — 144 particles are distributed onto 
233 periods of the substrate potential. 

The model can be enriched by considering the dynamics of the particles [TTl [12] . 
We include the kinetic energy of the particles by defining t — t/ro, where t is the time 
and tq = ^ m/ K . If we express all energies in units of Ka^, times in units of tq and 
lengths in units of Og, the Hamiltonian of the system becomes: 



where Ui is the position of particle i of the chain with the periodic boundary 
condition ^ ^ Ui. The resulting equations of motion can be written as a dynamical 
system of 2N first order ordinary differential equations: 





(2) 
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Figure 1. Top panel; a schematic overview of the Frenkel-Kontorova (FK) 
model. A chain of harmonically coupled particles is subjected to a periodic 
potential. Bottom panel: a sketch of the proposed phase diagram by Shinjo 
and Hirano (adapted from reference As the initial CM velocity is increased 

along the dashed line, the system supposedly first undergoes a transition from 
zero static friction (A < Ac) to kinetic friction and then to a super lubric state. 
Our results do not support this picture. 



— - = Ui+i + — 2ui — Asiii(27rui) . (3) 
dr 

We solve these equations numerically by using the fourth order Runge-Kutta (RK) 
method with a time step of to/150 and calculate also the position Q = (1/A^) X^i^i 
and velocity P — (l/N) X]i=i Vi of the CM. As the system is Hamiltonian (no heuristic 
damping) the total energy is conserved. Nevertheless, energy can be transferred from 
the centre of mass to the internal degrees of freedom, leading to the arrest of the chain 
in time. This effect can be interpreted as an effective friction. 

If there is no coupling between the chain and the substrate, A = 0, the dispersion 
relation of the phonon modes is the one of a harmonic chain, uik = 2|sin(fc/2)|, 
with LUk the frequency of the phonons and k — 2Tri/N the wavevectors in which 
i € {—N/2, N/2]. When all the particles move with the same velocity po, the atoms 
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slide over the periodic potential with the washboard frequency — 2ttpq. Strong 
resonances occur if ft is close to the frequency of the phonon mode with the wavevector 
q = 27rao/as or to its harmonics nq of multiples of the wavevector q. If A is small, these 
frequencies can be approximated by the frequencies aj„q of the phonons of the harmonic 
chain. The resonance phenomenon therefore occurs for CM velocities [13l [12] : 

npo — e N n > 1. (4) 

When these resonances occur, after an initial recurrent energy exchange pH I12| . 
the CM translational kinetic energy is eventually totally transformed into internal 
kinetic and potential energy of vibration in the chain. This energy transfer occurs 
because the CM motion with velocity P induces a modulation with the washboard 
frequency il that leads to parametric resonances with exponential growth of phonon 
modes with wavevector k whenever [T31 [T^l [T7] 

^ _ u}{k)+uj{mq- k) 
m 

for some integer m. From the above relation, one can derive a maximum value 
Pmax = ^max/i^Tr) = 4| cos(7rrg/2) | /(27r) 0.5254. Above this value, the transfer 
of energy from the CM to the phonons is blocked and the motion would be really 
superlubric. 

Lastly we note that, although the total energy is conserved, entropy increases 
when the energy is dispersed over the large number of internal degrees of freedom 
of the chain. Finite systems can show temporary decreases in the entropy, and, in 
principle, due to Poincare recurrence, the system could return close to the initial 
state with all energy in the CM momentum. In practice, however, for systems with 
a sufficiently large number of internal degrees of freedom, such a large decrease in 
entropy is highly improbable and occurs only on extremely long time scales. Indeed, 
we have never observed this reversal in the length of our simulations. 



3. Dissipation leading to thermal equilibrium 

In this section, we show the evolution of the CM velocity P for initial velocities close to 
the first four resonances cj„q and show how this process leads to thermal equilibrium. 
We take A — 0.05 « Ac/3 for which there is no static friction, namely we consider 
values of P on the vertical dashed line in figure [TJ For the initial positions of the 
particles we take a ground state of the system, calculated by minimizing the potential 
energy. We give all the particles the same initial velocity in order to have the CM 
velocity P equal to the velocity corresponding to the first four resonances (table [T]). 

In figure [2] we show the time evolution P{t) for the second resonance. We see 
an initial oscillatory phase, enlarged in the inset of figure [2Ia), where P periodically 
goes back to the initial value. After some time this recurrence mechanism breaks 
down and P decays with time. This effect, the onset of friction, was not reported in 
reference [IT], possibly, because of the use of a too small number of particles or too 
short integration times [12]. The mechanism for this behaviour, parametric excitation 
of acoustic phonons, has been analyzed in reference [T^]. A similar behaviour has also 
been found for sliding concentric nanotubes |16j . 

When the decay has finished, P jiggles around 0. From figure[5]Jc) we see that the 
CM kinetic energy has been transformed in equal amounts of kinetic and potential 
energy as expected classically for vibrations, a first indication that the system has 
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reached thermal equihbrium. For the other resonances, the evolution of P in time is 
qualitatively similar but with different time scales: the higher the order of resonance 
the longer the time before the start of the decay of P. For the fourth resonance, we had 
to use a slightly perturbed ground state to speed up this process. The increasingly long 
time it takes for the parametric growth of the phonon excitations to lead to friction 
makes it very difficult to study numerically the low velocity limit. We cannot rule 
out superlubric motion at very low velocities where the parametric growth of internal 
vibrations could occur either on time scales too long for the simulations or not occur 
at all. Notice that these results do not support the phase diagram [TT] reported in 
figure [TJ 

If, after the decay of P, the system has reached thermal equilibrium, equipartition 
should apply, namely the initial CM kinetic energy should be distributed over the N 
oscillators, each having /ce^eq- We use two independent methods to ascertain if this 
is the case. 

At thermal equilibrium, the probability density function of P{t) should be 
described by a Maxwell-Boltzmann (MB) distribution. In one dimension, the MB 
distribution coincides with a Gaussian distribution around 0, with variance — 
kBT^^/M, where M = Nm is the mass of the chain, fee is the Boltzmann constant 
and the equilibrium temperature. We calculate the probability density function 
from a normalized histogram of P, divided in 100 uniform bins, in the time interval in 
which P jiggles around 0. The result, shown for the second resonance in figure [^b) 
is indeed described well by a Gaussian with mean /z = 1.5 x 10~^ cls/tq and variance 
= 4.5 X 10"^ /tq. The estimated values for the equilibrium temperature are 
given in table [TJ 

Another estimate for the temperature can be given by the amplitudes, Xk, of the 
phonon modes [T7]: J2k=i ^^^ki^k) = k^^T^q"^^, where the brackets denote a time 
average and TJjJ^°'^ is the equilibrium temperature corresponding to the distribution 
of the phonon amplitudes. In figure [21(d) we show, for the second resonance, that the 
amplitudes of the phonon modes are proportional to 1 /o;^ as expected at equilibrium. 

The behaviour of the distribution of P and the phonon modes is qualitatively the 
same for the other resonances. In table [T] we compare the temperatures Tcq with the 
values obtained from the distribution of P, T^^, and from the phonon amplitudes, 
TP|^°". In order to define high and low temperature regimes, all temperatures are 
given in units of the value of A used in this work, 0.05. In this way we can distinguish 
high (T « 1 A) and low temperature (T << 1 A) regimes. In the next section, we also 
study the generic value T — 0.16 A, between first and second resonance, for which the 
distribution of P was checked and found to be a Gaussian. The agreement between 
the different estimates demonstrates that the system has reached thermal equilibrium 
in the studied range of velocities around the first four resonances. 

4. Emergence of viscous damping 

After the system has reached thermal equilibrium, as described in the previous section, 
we give the CM an initial velocity P{t = 0) and study how it evolves with time. If the 
friction is purely viscous, Fhic — —rjMP, as is often assumed in the study of frictional 
dynamics |13l 1181 119] , the sliding velocity decays exponentially with time with a time 
scale determined by the friction coefficient 77. This type of viscous behaviour can be 
expected on general theoretical grounds for many-particle systems |20[ 121] . 
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Table 1. For each resonance n we give the initial velocity per particle pQ. t- 
range gives the time range in which the probability density function of P and 
the averages of the phonon amplitudes have been calculated. The absolute value 
of the averages (fi) of the Gaussian fits of the histograms of P is given in the 
fourth column in 10~*as/ro. They are a measure of the deviation of the system 
from a perfect thermal equilibrium. The last three columns are the temperatures 
obtained from the initial kinetic energy (Tcq), the distribution of P (T^^) and 

the phonon modes (TjjJ^""). 



n 




T-range (ro) 




Tcq (A) 


T^-(A) 


TPf°"(A) 


1 


0.2966 


(5 - 405) X 10-* 


6.5 ±2 


0.88 


0.89 ±0.01 


0.91 ±0.04 


2 


0.1075 


(25-425) X 10"* 


1.5±0.7 


0.12 


0.13 ±0.01 


0.13 ±0.01 


3 


0.04692 


(4- 16) X 10^ 


3.1 ± 1 


0.022 


0.026 ±0.002 


0.030 ±0.004 


4 


0.07927 


(4- 8) X 10^ 


2.0±0.7 


0.063 


0.055 ±0.002 


0.074 ±0.004 



0.12 
0.1 
0.08 
0.06 
0.04 
0.02 


-0.02 



(a) 






m 


















100 200 300 




2e+04 4e+04 6e+04 8e+04 10e+04 



-0.02 -0.01 0.01 0.02 

P(as/%) 




26-1-04 4e-i-04 66-1-04 86-1-04 106-1-04 




Figure 2. Time evolution of the system for the second resonance, (a) P as 
a function of r in the region in which P decays with time. The inset shows the 
initial recurrent mechanism, (b) Probability density function of P from data of 
the P-T evolution (points) and the Gaussian fit of this function (line), (c) E as 
a function of r for the second order of resonance. The total kinetic energy and 
the increase in potential energy relative to the minimum configuration are shown. 
The dashed vertical line (red) is the time at which equipartition is reached, (d) 
The amplitudes of the phonon modes < x| > as a function of mode number k 
for the second order of resonance. The data points are plotted in green and the 
fit according to l/t^^^ is in blue. The deviation of the behaviour around 

fc = 55 = —q is due to the approximated incommensurability of the system. For 
the time ranges used for averaging in (b) and (d), see Table [T] 



We start from initial conditions obtained by choosing many (typically hundred) 
configurations of the system in thermal equilibrium at temperature Tcq in the r- 
range given in table [T] We increase the CM velocity P by adding a uniform initial 
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Figure 3. Illustration of the way to calculate r]. The evolution of {P{t)) with 
respect to time after the addition of a velocity p+ = 0.09 to all the particles 
for T = 0.88 A, averaged over 500 initial configurations chosen with a spacing 
in time of 2000 tq from the trajectory starting from the first resonance, with 
T = 0.88 A. The fit is an exponential function, (P(r))(as/To) = 0.0923 exp(-7jr) 
with 7? = 0.00194(1/to). 



sliding velocity p"*" to all the particles. The systems are then evolved with time and 
the time dependence of the averaged CM velocity {P{t)) is obtained by averaging 
over all initial configurations. This averaging is necessary, as, in such small systems, 
fluctuations are large. In figure[3]we show (-P(r)) for p+ = 0.09 Os/tq, averaged over 
500 configurations at T = 0.88 A. The average CM velocity decays exponentially, 
{P{T)){ajTo) = 0.0923 exp(-7/r) with a friction coefficient = 1.94 x 10-3(1/to). 

In this system, it is the microscopic dynamics that give rise to the friction 
and determine the value of the friction coefficient, without the need for imposing 
a phenomenological friction term. A related result for damping of phonon modes in 
relation to thermostats in crystalline materials was obtained and tested in the context 
of friction [22] . 

The time range used for the calculation of rj needs to be not too short because 
of the reliability of the fit, nor too long because of the change of temperature of the 
system which occurs during the time range. We choose for the high temperature case, 
the time range r = — 500 tq and for the low temperature case, the time range 
T = — 1000 To, to determine r] and its error bar. Therefore, we do not consider the 
effect of different time ranges for the fit, although we have checked that the behaviour 
of rj remains qualitatively the same. 

The behaviour of 77 as a function of sliding velocity is shown in figure S) For 
all T, in the limit ^ 0, rj becomes constant (as indicated by the arrows in 
figure [4]), implying that the viscous friction, Ffj-ic, goes to zero and superlubricity is 
conserved. The same conclusion has been drawn in |23| for the Tomlinson model. The 
friction coefficient 77 grows with temperature. This is not surprising, because thermal 
fluctuations distort the incommensurate contact between the surfaces. It is known that 
friction decreases with temperature, particularly in the stick slip regime [241 123] where 
the 'slips' from one minimum of the potential energy to the other can be activated 
by thermal energy. The attempt frequency of these processes is determined by the 
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damping rj that is usually taken to be temperature independent. 

The vertical dashed lines in figure H] are the velocities corresponding to the 
first four resonances of the system (equation 0). For the high temperature case 
in figure Hta), the friction coefficient has a peak at frequencies just above the first 
resonance, ujq and no pronounced features at the lower three resonances. For the low 
temperature case, figure |4ljb), the influence of the lower resonances on the friction 
characteristics becomes apparent. Again the peaks in 77 are shifted to the right. In 
addition, the influence of temperature on the broadening of the peaks is visible, the 
higher the temperature, the broader the peak. For the lowest studied resonance, the 
fourth, there is no clear peak of rj at any of the studied temperatures, whereas for 
the third resonance it can be seen at the lowest studied temperature. Nevertheless, at 
sliding velocities near the resonances, superlubricity is suppressed for all temperatures. 

In section [5] we have shown that at T = there is a maximum resonance velocity, 
Pmax ~ 0.5254. For initial velocities sufficiently far above this value, there can be no 
effective damping at T = 0. At finite temperature, thermal fluctuations distort the 
incommensurate geometry, and nonzero effective damping can still exist even for high 
velocities. We have indeed observed a slow decay of velocity in our simulations at 
p+ = 0.75 and T = 0.88A. 

Our results can be compared to the velocity-force characteristics calculated by 
Strunz and Elmer |13) for a driven and damped FK- model at T = 0. The inherent 
friction in their system, calculated by subtracting the imposed damping term from 
the total force, shows an effect of the resonances similar to our results of figure ID A 
difference between our results and those of Strunz and Elmer is the location of the 
peaks in the friction characteristics. In reference [13] these peaks are situated at the 
values of the resonances calculated from equation Q , while in our situation the peaks 
are shifted towards higher values of the sliding velocity. A possible explanation for 
this shift towards the right in our undamped case, is the influence of the use of the 
approximation A = for the resonances equation ([4]). Indeed, for smaller coupling, 
A = 0.005 and T = 8.8 A, we have observed a smaller shift for the first resonance. 

Also Strunz and Elmer used this approximation, but they studied the so called 
uniform sliding state in a damped system, whereas we solve the full dynamics. Our 
results confirm the infiuence of the resonances on the friction found in reference |13j . 
and show that this aspect of the model is not limited to T = nor to uniform sliding 
states of a damped driven system. 

5. Estimate of parameters 

Our ID model is highly simplified and focuses on the geometrical incommensurability 
which is essential for superlubricity. In real systems also other aspects like 
roughness [6] and electronic friction [7j [25l [26] may play a role. Nevertheless, 
we attempt to find possible experimental situations in which our results could be 
confirmed. An example of a well studied system with incommensurate contacting 
surfaces is a layer of Xe atoms sliding on an Ag(lll) substrate [261 [27l [28]. The 
ratio between the lattice parameters of the Xe monolayer and the substrate is 
ao/tts = (4.55 X 10-i°m/2.892 x lO'^^m) = 1.57 « Tg. The mass of the Xe atoms is 
m = 2.16 X 10^^^ kg. The spring constant between the Xe atoms can be obtained 
from a Taylor expansion of the Lennard- Jones potential with parameters for Xe-Xe 
around its minimum [26], K = 1.103 J/m^. Consequently, the typical scale of the 
velocity is Qs/to = 6.5 x 10^ m/s. In a similar way Uq is estimated to be in the range 
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Figure 4. 77 as a function of p+ for different temperatures, (a) the high 
temperature case, T = 0.88 A, and (b) the low temperature cases. The vertical 
dashed lines are the values of p+ corresponding to the first four resonances 
(equation (|4]l). The peak resulting from the resonance broadens at higher 
temperatures. Each friction coefficient was obtained from averaging over at least 
too initial configurations separated 2000 tq in time, tj has been obtained from 
T = - 500 To in (a) and from r = - fOOO tq in (b). 



1 - 5 X 10-21 J [261 [29], which gives 0.01 < A < 0.05 < Ac close to the value of A we 
used in the simulations. 

The velocity where the effective friction becomes a strongly nonlinear function 
of velocity is rather high, order of magnitudes too high for AFM and at the limit 
of reachable velocities for QCM. Therefore it is important to realize that in the low 
velocity limit, friction grows smoothly with velocity from the vanishing static friction, 
contrary to the Shinjo-Hirano model of figure [Tj It would be interesting however to 
explore also the transition from low viscous friction to high friction in correspondence 
of the resonances described before. The friction characteristics for the third resonance 
could possibly be seen at sliding velocities of approximately 30 m/s, for temperature 
in the range of about 8 K for Xe sliding on Ag(lll). The second resonance, at around 
70 m/s, would be visible up to higher temperatures, at least around 50 K. 
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Another interesting experimental setting for the FK model is in cohoidal 
suspensions, where atomic-scale systems can be reproduced on larger scales, with 
colloidal particles playing the role of atoms. Recently, the static FK model was studied 
experimentally in such systems [30 . However, as colloidal crystals have over-damped 
decay of phonon modes [31] , the present results cannot easily be scaled up in the same 
way. 

6. Summary 

We have numerically calculated effective friction coefficients as a function of 
temperature and velocity in the FK-model without any driving forces or 
phenomenological damping terms. We have critically reviewed previous findings 
and shown how initial kinetic energy of the centre of mass is converted into heat 
by excitation of the internal motion of the chain, bringing the system to thermal 
equilibrium. We have then studied the temperature and velocity dependence of 
friction and found that it can be described as an effective viscous damping. The 
friction coefficient however, shows a peculiar velocity dependence with peaks at sliding 
velocities close to resonances between the phonons of the chain and wavevectors 
related to the modulation induced by the incommensurate periodic potential. We 
argue that these peaks in the friction coefficient could be observed experimentally. 
Moreover, we find that for incommensurate contacts, the effective damping increases 
with temperature. 

Despite the fact that both sliding velocity and temperature can suppress 
superlubricity in incommensurate systems, the fact that the friction coefficient remains 
finite at low velocities, makes that superlubricity is preserved in the important limit 
of very low velocity probed by AFM experiments. 
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